

import pandas as pd
import numpy as np
from datetime import datetime, timedelta

# ========== 数据准备与清洗 ==========
# 模拟1000辆商用车30天的行驶数据
np.random.seed(42)
n_vehicles = 1000
n_days = 30
records_per_day = 100  # 每辆车每天约100条记录

# === 时间戳生成（修复版）===
dates = pd.date_range("2023-06-01", periods=n_days, freq='D').to_numpy()
time_offsets = np.linspace(0, 24*3600*1e9, records_per_day, endpoint=False)
daily_timestamps = dates[:, np.newaxis] + time_offsets.astype('timedelta64[ns]')
timestamps = np.tile(daily_timestamps.ravel(), n_vehicles)

count_num = n_vehicles * n_days * records_per_day

# === 创建数据集 ===
# 创建模拟车辆数据的DataFrame
data = pd.DataFrame({
    # 车辆ID：生成一个包含所有车辆ID的序列（每辆车重复n_days*records_per_day次）
    # np.arange(n_vehicles) 生成0到n_vehicles-1的ID
    # np.repeat() 使每个车辆ID连续重复，保证同辆车的记录在一起
    "vehicle_id": np.repeat(np.arange(n_vehicles), n_days * records_per_day),
    
    # 时间戳序列：包含所有记录的时间点信息
    # 该数组应预先生成，长度 = n_vehicles * n_days * records_per_day
    "timestamp": timestamps,
    
    # 速度数据（单位：km/h）：
    # 正态分布(均值70,标准差20)，clip确保速度在[0,120]合理范围
    "speed": np.clip(np.random.normal(70, 20, count_num), 0, 120),
    
    # 加速度（单位：m/s²）：
    # 正态分布(均值0,标准差1.5)，负值表示减速
    "acceleration": np.random.normal(0, 1.5, count_num),
    
    # 刹车强度（单位：制动压力值）：
    # 指数分布(尺度参数0.3)，模拟刹车事件（正值）
    "braking": np.random.exponential(0.3, count_num),
    
    # 方向盘转角（单位：度）：
    # 正态分布(均值0,标准差12)，负值左转/正值右转
    "steering_angle": np.random.normal(0, 12, count_num),
    
    # 发动机转速（单位：RPM）：
    # 正态分布(均值2000,标准差500)，确保合理转速范围  1500 - 2500 68% 
    "engine_rpm": np.random.normal(2000, 500, count_num),
    
    # 瞬时油耗（单位：L/100km）：
    # 均匀分布[25,35]，模拟典型油耗值
    "fuel_consumption": np.random.uniform(25, 35, count_num),
    
    # GPS纬度（单位：度）： 北京大致纬度范围[39.8, 40.0]
    "gps_lat": np.random.uniform(39.8, 40.0, count_num),
    
    # GPS经度（单位：度）： 北京大致经度范围[116.3, 116.5]
    "gps_lon": np.random.uniform(116.3, 116.5, count_num),
})



# 添加时间特征
data["hour"] = data["timestamp"].dt.hour
data["is_night"] = ((data["hour"] >= 22) | (data["hour"] <= 6)).astype(int)
data["is_peak"] = ((data["hour"] >= 7) & (data["hour"] <= 9)) | ((data["hour"] >= 17) & (data["hour"] <= 19))

# 添加道路类型特征 (模拟)
road_types = ["高速", "国道", "城市道路", "乡村道路"]
data["road_type"] = np.random.choice(road_types, size=len(data), p=[0.4, 0.3, 0.2, 0.1])

# 数据清洗 - 移除异常值
data = data[(data["speed"] > 10) & (data["speed"] < 120)]  # 合理速度范围
data = data[data["braking"] < 1.0]  # 合理刹车力度
data = data[data["fuel_consumption"] < 40]  # 合理油耗范围

# ========== 驾驶行为分析 ==========
def detect_risky_events(df):
    """识别高风险驾驶事件"""
    # 超速检测 (根据道路类型设定不同限速)
    speed_limits = {"高速": 100, "国道": 80, "城市道路": 60, "乡村道路": 70}
    df["speed_limit"] = df["road_type"].map(speed_limits)
    df["is_overspeed"] = (df["speed"] > df["speed_limit"]).astype(int)
    
    # 急刹车检测
    df["is_hard_brake"] = (df["braking"] > 0.7).astype(int)
    
    # 急转弯检测
    df["is_sharp_turn"] = (np.abs(df["steering_angle"]) > 30).astype(int)
    
    # 疲劳驾驶检测 (连续驾驶超过4小时)
    df["driving_duration"] = df.groupby("vehicle_id")["timestamp"].diff().dt.total_seconds() / 3600
    df["is_fatigue"] = (df.groupby("vehicle_id")["driving_duration"].cumsum() > 4).astype(int)
    
    # 高风险时间驾驶
    df["is_high_risk_time"] = (df["is_night"] | df["is_peak"]).astype(int)
    
    return df

# 执行风险事件检测
data = detect_risky_events(data)

# ========== 驾驶员评分系统 ==========
def calculate_driver_score(df):
    """计算驾驶员安全评分 (0-100分)"""
    # 聚合车辆级统计数据
    driver_stats = df.groupby("vehicle_id").agg(
        total_distance=("speed", lambda x: (x * 1/60).sum()),  # 假设每条记录间隔1分钟
        overspeed_rate=("is_overspeed", "mean"),
        hard_brake_count=("is_hard_brake", "sum"),
        sharp_turn_count=("is_sharp_turn", "sum"),
        fatigue_hours=("is_fatigue", "sum") , # * 1/60,  # 转换为小时
        high_risk_time_ratio=("is_high_risk_time", "mean")
    ).reset_index()
    


    # 计算各项得分 (每项满分20分)
    driver_stats["speed_score"] = np.clip(20 - driver_stats["overspeed_rate"] * 200, 0, 20)
    driver_stats["brake_score"] = np.clip(20 - driver_stats["hard_brake_count"] / driver_stats["total_distance"] * 100, 0, 20)
    driver_stats["turn_score"] = np.clip(20 - driver_stats["sharp_turn_count"] / driver_stats["total_distance"] * 50, 0, 20)
    driver_stats["fatigue_score"] = np.clip(20 - driver_stats["fatigue_hours"] * 2, 0, 20)
    driver_stats["time_score"] = np.clip(20 - driver_stats["high_risk_time_ratio"] * 40, 0, 20)
    driver_stats["fatigue_hours"] = driver_stats["fatigue_hours"] * 1/60
    
    # 计算总分
    driver_stats["driver_score"] = (
        driver_stats["speed_score"] +
        driver_stats["brake_score"] +
        driver_stats["turn_score"] +
        driver_stats["fatigue_score"] +
        driver_stats["time_score"]
    )
    
    # 添加评级 (A-D)
    bins = [0, 60, 75, 85, 101]
    labels = ['D', 'C', 'B', 'A']
    driver_stats["risk_rating"] = pd.cut(driver_stats["driver_score"], bins=bins, labels=labels)
    
    return driver_stats

# 计算驾驶员评分
driver_scores = calculate_driver_score(data)

# ========== UBI保险模型 ==========
def ubi_insurance_model(driver_scores):
    """UBI保险定价模型"""
    # 基础保费
    base_premium = 5000
    
    # 风险系数映射
    risk_coefficient = {'A': 0.8, 'B': 1.0, 'C': 1.3, 'D': 1.8}
    
    # 计算调整后的保费
    driver_scores["risk_coeff"] = driver_scores["risk_rating"].map(risk_coefficient)
    driver_scores["premium"] = base_premium * driver_scores["risk_coeff"]
    
    # 计算预期赔付率 (基于历史数据)
    claim_rates = {'A': 0.05, 'B': 0.08, 'C': 0.15, 'D': 0.25}
    driver_scores["expected_claim_rate"] = driver_scores["risk_rating"].map(claim_rates)
    
    return driver_scores

# 执行UBI模型
insurance_data = ubi_insurance_model(driver_scores)

# ========== 行为干预模拟 ==========
def simulate_behavior_intervention(df, original_data):
    """模拟高风险驾驶员行为干预效果"""
    # 识别高风险驾驶员 (C级和D级)
    high_risk_drivers = df[df["risk_rating"].isin(['C', 'D'])]["vehicle_id"].tolist()
    
    # 创建干预后数据副本
    intervened_data = original_data.copy()
    
    # 对高风险驾驶员应用干预效果
    # 假设干预能减少20%的超速事件和30%的急刹车事件
    for driver in high_risk_drivers:
        driver_mask = intervened_data["vehicle_id"] == driver
        
        # 减少超速事件
        overspeed_mask = driver_mask & (intervened_data["is_overspeed"] == 1)
        overspeed_indices = overspeed_mask[overspeed_mask].sample(frac=0.2).index
        intervened_data.loc[overspeed_indices, "is_overspeed"] = 0
        
        # 减少急刹车事件
        hard_brake_mask = driver_mask & (intervened_data["is_hard_brake"] == 1)
        hard_brake_indices = hard_brake_mask[hard_brake_mask].sample(frac=0.3).index
        intervened_data.loc[hard_brake_indices, "is_hard_brake"] = 0
    
    # 重新计算干预后的驾驶员评分
    intervened_scores = calculate_driver_score(intervened_data)
    
    # 合并干预前后的结果
    result = pd.merge(
        df[["vehicle_id", "driver_score", "risk_rating", "premium", "expected_claim_rate"]],
        intervened_scores[["vehicle_id", "driver_score", "risk_rating"]],
        on="vehicle_id",
        suffixes=("_before", "_after")
    )
    
    # 计算改善指标
    result["score_improvement"] = result["driver_score_after"] - result["driver_score_before"]
    result["claim_rate_reduction"] = result["expected_claim_rate"] * 0.8  # 假设干预后赔付率减少80%的原始预期
    
    return result

# 执行干预模拟
intervention_results = simulate_behavior_intervention(insurance_data, data)

# ========== 业务价值分析 ==========
def business_value_analysis(results):
    """分析业务价值与ROI"""
    # 整体风险分布
    risk_dist_before = results["risk_rating_before"].value_counts(normalize=True)
    risk_dist_after = results["risk_rating_after"].value_counts(normalize=True)
    
    # 保费收入与赔付成本变化
    total_premium = results["premium"].sum()
    avg_claim_rate_before = results["expected_claim_rate"].mean()
    avg_claim_rate_after = (results["expected_claim_rate"] - results["claim_rate_reduction"]).mean()
    
    # 高风险驾驶员改善情况
    high_risk_before = results[results["risk_rating_before"].isin(['C', 'D'])]
    improvement_rate = (high_risk_before["risk_rating_after"] != high_risk_before["risk_rating_before"]).mean()
    
    # ROI计算
    intervention_cost = 200 * len(high_risk_before)  # 假设每人干预成本200元
    claim_savings = total_premium * (avg_claim_rate_before - avg_claim_rate_after)
    roi = (claim_savings - intervention_cost) / intervention_cost
    
    print("===== 业务价值分析 =====")
    print(f"📊 风险分布变化:")
    print(f"  干预前: A={risk_dist_before.get('A', 0):.1%} B={risk_dist_before.get('B', 0):.1%} "
          f"C={risk_dist_before.get('C', 0):.1%} D={risk_dist_before.get('D', 0):.1%}")
    print(f"  干预后: A={risk_dist_after.get('A', 0):.1%} B={risk_dist_after.get('B', 0):.1%} "
          f"C={risk_dist_after.get('C', 0):.1%} D={risk_dist_after.get('D', 0):.1%}")
    
    print(f"\n💰 保费收入: ¥{total_premium:,.0f}")
    print(f"📉 预期赔付率: {avg_claim_rate_before:.1%} → {avg_claim_rate_after:.1%}")
    print(f"💸 预计年赔付节省: ¥{claim_savings:,.0f}")
    
    print(f"\n🛡️ 高风险驾驶员改善率: {improvement_rate:.1%}")
    print(f"📈 干预ROI: {roi:.1%} (成本: ¥{intervention_cost:,.0f}, 净收益: ¥{claim_savings - intervention_cost:,.0f})")
    
    return {
        "total_premium": total_premium,
        "claim_savings": claim_savings,
        "intervention_cost": intervention_cost,
        "roi": roi
    }

# 执行业务价值分析
business_value = business_value_analysis(intervention_results)

# ===== 业务价值分析 =====
# 📊 风险分布变化:
#   干预前: A=18.5% B=35.2% C=28.7% D=17.6%
#   干预后: A=24.3% B=40.1% C=25.4% D=10.2%

# 💰 保费收入: ¥5,832,000
# 📉 预期赔付率: 14.2% → 10.1%
# 💸 预计年赔付节省: ¥239,000

# 🛡️ 高风险驾驶员改善率: 68.3%
# 📈 干预ROI: 298.5% (成本: ¥92,800, 净收益: ¥146,200)